#include <cctype>
#include <cmath>

#include "../gmx_blas.h"
#include "../gmx_lapack.h"

#include "gromacs/utility/real.h"

void F77_FUNC(sbdsqr, SBDSQR)(const char* uplo,
                              int*        n,
                              int*        ncvt,
                              int*        nru,
                              int*        ncc,
                              float*      d__,
                              float*      e,
                              float*      vt,
                              int*        ldvt,
                              float*      u,
                              int*        ldu,
                              float*      c__,
                              int*        ldc,
                              float*      work,
                              int*        info)
{
    const char xuplo = std::toupper(*uplo);
    int        c_dim1, c_offset, u_dim1, u_offset, vt_dim1, vt_offset, i__1, i__2;
    float      r__1, r__2, r__3, r__4;
    float      c_b15 = -.125;

    int   c__1  = 1;
    float c_b49 = 1.f;
    float c_b72 = -1.f;

    float f, g, h__;
    int   i__, j, m;
    float r__, cs;
    int   ll;
    float sn, mu;
    int   nm1, nm12, nm13, lll;
    float eps, sll, tol, abse;
    int   idir;
    float abss;
    int   oldm;
    float cosl;
    int   isub, iter;
    float unfl, sinl, cosr, smin, smax, sinr;
    float oldcs;
    int   oldll;
    float shift, sigmn, oldsn = 0.;
    int   maxit;
    float sminl;
    float sigmx;
    int   lower;
    float sminoa;
    float thresh;
    int   rotate;
    float tolmul;
    int   itmp1, itmp2;

    --d__;
    --e;
    vt_dim1   = *ldvt;
    vt_offset = 1 + vt_dim1;
    vt -= vt_offset;
    u_dim1   = *ldu;
    u_offset = 1 + u_dim1;
    u -= u_offset;
    c_dim1   = *ldc;
    c_offset = 1 + c_dim1;
    c__ -= c_offset;
    --work;

    *info = 0;

    itmp1 = (*n > 1) ? *n : 1;
    itmp2 = (*nru > 1) ? *nru : 1;

    lower = (xuplo == 'L');
    if ((xuplo != 'U') && !lower)
    {
        *info = -1;
    }
    else if (*n < 0)
    {
        *info = -2;
    }
    else if (*ncvt < 0)
    {
        *info = -3;
    }
    else if (*nru < 0)
    {
        *info = -4;
    }
    else if (*ncc < 0)
    {
        *info = -5;
    }
    else if (((*ncvt == 0) && (*ldvt < 1)) || ((*ncvt > 0) && (*ldvt < itmp1)))
    {
        *info = -9;
    }
    else if (*ldu < itmp2)
    {
        *info = -11;
    }
    else if (((*ncc == 0) && (*ldc < 1)) || ((*ncc > 0) && (*ldc < itmp1)))
    {
        *info = -13;
    }
    if (*info != 0)
    {
        return;
    }
    if (*n == 0)
    {
        return;
    }
    if (*n == 1)
    {
        goto L160;
    }

    rotate = *ncvt > 0 || *nru > 0 || *ncc > 0;

    if (!rotate)
    {
        F77_FUNC(slasq1, SLASQ1)(n, &d__[1], &e[1], &work[1], info);
        return;
    }

    nm1  = *n - 1;
    nm12 = nm1 + nm1;
    nm13 = nm12 + nm1;
    idir = 0;

    eps  = GMX_FLOAT_EPS;
    unfl = GMX_FLOAT_MIN / GMX_FLOAT_EPS;

    if (lower)
    {
        i__1 = *n - 1;
        for (i__ = 1; i__ <= i__1; ++i__)
        {
            F77_FUNC(slartg, SLARTG)(&d__[i__], &e[i__], &cs, &sn, &r__);
            d__[i__]        = r__;
            e[i__]          = sn * d__[i__ + 1];
            d__[i__ + 1]    = cs * d__[i__ + 1];
            work[i__]       = cs;
            work[nm1 + i__] = sn;
        }

        if (*nru > 0)
        {
            F77_FUNC(slasr, SLASR)("R", "V", "F", nru, n, &work[1], &work[*n], &u[u_offset], ldu);
        }
        if (*ncc > 0)
        {
            F77_FUNC(slasr, SLASR)("L", "V", "F", n, ncc, &work[1], &work[*n], &c__[c_offset], ldc);
        }
    }

    r__3 = 100.f, r__4 = std::pow(static_cast<float>(GMX_FLOAT_EPS), c_b15);
    r__1 = 10.f, r__2 = (r__3 < r__4) ? r__3 : r__4;
    tolmul = (r__1 > r__2) ? r__1 : r__2;
    tol    = tolmul * eps;
    smax   = 0.f;
    i__1   = *n;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        r__2 = smax, r__3 = (r__1 = d__[i__], std::abs(r__1));
        smax = (r__2 > r__3) ? r__2 : r__3;
    }
    i__1 = *n - 1;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        r__2 = smax, r__3 = (r__1 = e[i__], std::abs(r__1));
        smax = (r__2 > r__3) ? r__2 : r__3;
    }
    sminl = 0.f;
    if (tol >= 0.f)
    {
        sminoa = std::abs(d__[1]);
        if (sminoa == 0.f)
        {
            goto L50;
        }
        mu   = sminoa;
        i__1 = *n;
        for (i__ = 2; i__ <= i__1; ++i__)
        {
            mu = (r__2 = d__[i__], std::abs(r__2)) * (mu / (mu + (r__1 = e[i__ - 1], std::abs(r__1))));
            sminoa = (sminoa < mu) ? sminoa : mu;
            if (sminoa == 0.f)
            {
                goto L50;
            }
        }
    L50:
        sminoa /= std::sqrt((float)(*n));
        r__1 = tol * sminoa, r__2 = *n * 6 * *n * unfl;
        thresh = (r__1 > r__2) ? r__1 : r__2;
    }
    else
    {
        r__1 = std::abs(tol) * smax, r__2 = *n * 6 * *n * unfl;
        thresh = (r__1 > r__2) ? r__1 : r__2;
    }
    maxit = *n * 6 * *n;
    iter  = 0;
    oldll = -1;
    oldm  = -1;
    m     = *n;

L60:

    if (m <= 1)
    {
        goto L160;
    }
    if (iter > maxit)
    {
        goto L200;
    }

    if (tol < 0.f && (r__1 = d__[m], std::abs(r__1)) <= thresh)
    {
        d__[m] = 0.f;
    }
    smax = (r__1 = d__[m], std::abs(r__1));
    smin = smax;
    i__1 = m - 1;
    for (lll = 1; lll <= i__1; ++lll)
    {
        ll   = m - lll;
        abss = (r__1 = d__[ll], std::abs(r__1));
        abse = (r__1 = e[ll], std::abs(r__1));
        if (tol < 0.f && abss <= thresh)
        {
            d__[ll] = 0.f;
        }
        if (abse <= thresh)
        {
            goto L80;
        }
        smin = (smin < abss) ? smin : abss;
        r__1 = (smax > abss) ? smax : abss;
        smax = (r__1 > abse) ? r__1 : abse;
    }
    ll = 0;
    goto L90;
L80:
    e[ll] = 0.f;
    if (ll == m - 1)
    {
        --m;
        goto L60;
    }
L90:
    ++ll;
    if (ll == m - 1)
    {
        F77_FUNC(slasv2, SLASV2)
        (&d__[m - 1], &e[m - 1], &d__[m], &sigmn, &sigmx, &sinr, &cosr, &sinl, &cosl);
        d__[m - 1] = sigmx;
        e[m - 1]   = 0.f;
        d__[m]     = sigmn;
        if (*ncvt > 0)
        {
            F77_FUNC(srot, SROT)
            (ncvt, &vt[m - 1 + vt_dim1], ldvt, &vt[m + vt_dim1], ldvt, &cosr, &sinr);
        }
        if (*nru > 0)
        {
            F77_FUNC(srot, SROT)
            (nru, &u[(m - 1) * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], &c__1, &cosl, &sinl);
        }
        if (*ncc > 0)
        {
            F77_FUNC(srot, SROT)
            (ncc, &c__[m - 1 + c_dim1], ldc, &c__[m + c_dim1], ldc, &cosl, &sinl);
        }
        m += -2;
        goto L60;
    }
    if (ll > oldm || m < oldll)
    {
        if ((r__1 = d__[ll], std::abs(r__1)) >= (r__2 = d__[m], std::abs(r__2)))
        {
            idir = 1;
        }
        else
        {
            idir = 2;
        }
    }
    if (idir == 1)
    {

        if ((std::abs(e[m - 1]) <= std::abs(tol) * std::abs(d__[m]))
            || (tol < 0.0 && std::abs(e[m - 1]) <= thresh))
        {
            e[m - 1] = 0.f;
            goto L60;
        }
        if (tol >= 0.f)
        {
            mu    = (r__1 = d__[ll], std::abs(r__1));
            sminl = mu;
            i__1  = m - 1;
            for (lll = ll; lll <= i__1; ++lll)
            {
                if ((r__1 = e[lll], std::abs(r__1)) <= tol * mu)
                {
                    e[lll] = 0.f;
                    goto L60;
                }
                mu = (r__2 = d__[lll + 1], std::abs(r__2))
                     * (mu / (mu + (r__1 = e[lll], std::abs(r__1))));
                sminl = (sminl < mu) ? sminl : mu;
            }
        }
    }
    else
    {
        if ((std::abs(e[ll]) <= std::abs(tol) * std::abs(d__[ll]))
            || (tol < 0.0 && std::abs(e[ll]) <= thresh))
        {
            e[ll] = 0.f;
            goto L60;
        }
        if (tol >= 0.f)
        {
            mu    = (r__1 = d__[m], std::abs(r__1));
            sminl = mu;
            i__1  = ll;
            for (lll = m - 1; lll >= i__1; --lll)
            {
                if ((r__1 = e[lll], std::abs(r__1)) <= tol * mu)
                {
                    e[lll] = 0.f;
                    goto L60;
                }
                mu = (r__2 = d__[lll], std::abs(r__2)) * (mu / (mu + (r__1 = e[lll], std::abs(r__1))));
                sminl = (sminl < mu) ? sminl : mu;
            }
        }
    }
    oldll = ll;
    oldm  = m;

    r__1 = eps, r__2 = tol * .01f;
    if (tol >= 0.f && *n * tol * (sminl / smax) <= ((r__1 > r__2) ? r__1 : r__2))
    {
        shift = 0.f;
    }
    else
    {
        if (idir == 1)
        {
            sll = (r__1 = d__[ll], std::abs(r__1));
            F77_FUNC(slas2, SLAS2)(&d__[m - 1], &e[m - 1], &d__[m], &shift, &r__);
        }
        else
        {
            sll = (r__1 = d__[m], std::abs(r__1));
            F77_FUNC(slas2, SLAS2)(&d__[ll], &e[ll], &d__[ll + 1], &shift, &r__);
        }
        if (sll > 0.f)
        {
            r__1 = shift / sll;
            if (r__1 * r__1 < eps)
            {
                shift = 0.f;
            }
        }
    }
    iter = iter + m - ll;
    if (shift == 0.f)
    {
        if (idir == 1)
        {
            cs    = 1.f;
            oldcs = 1.f;
            i__1  = m - 1;
            for (i__ = ll; i__ <= i__1; ++i__)
            {
                r__1 = d__[i__] * cs;
                F77_FUNC(slartg, SLARTG)(&r__1, &e[i__], &cs, &sn, &r__);
                if (i__ > ll)
                {
                    e[i__ - 1] = oldsn * r__;
                }
                r__1 = oldcs * r__;
                r__2 = d__[i__ + 1] * sn;
                F77_FUNC(slartg, SLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
                work[i__ - ll + 1]        = cs;
                work[i__ - ll + 1 + nm1]  = sn;
                work[i__ - ll + 1 + nm12] = oldcs;
                work[i__ - ll + 1 + nm13] = oldsn;
            }
            h__      = d__[m] * cs;
            d__[m]   = h__ * oldcs;
            e[m - 1] = h__ * oldsn;
            if (*ncvt > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[ll + vt_dim1], ldvt);
            }
            if (*nru > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 + 1], &u[ll * u_dim1 + 1], ldu);
            }
            if (*ncc > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 + 1], &c__[ll + c_dim1], ldc);
            }
            if ((r__1 = e[m - 1], std::abs(r__1)) <= thresh)
            {
                e[m - 1] = 0.f;
            }
        }
        else
        {
            cs    = 1.f;
            oldcs = 1.f;
            i__1  = ll + 1;
            for (i__ = m; i__ >= i__1; --i__)
            {
                r__1 = d__[i__] * cs;
                F77_FUNC(slartg, SLARTG)(&r__1, &e[i__ - 1], &cs, &sn, &r__);
                if (i__ < m)
                {
                    e[i__] = oldsn * r__;
                }
                r__1 = oldcs * r__;
                r__2 = d__[i__ - 1] * sn;
                F77_FUNC(slartg, SLARTG)(&r__1, &r__2, &oldcs, &oldsn, &d__[i__]);
                work[i__ - ll]        = cs;
                work[i__ - ll + nm1]  = -sn;
                work[i__ - ll + nm12] = oldcs;
                work[i__ - ll + nm13] = -oldsn;
            }
            h__     = d__[ll] * cs;
            d__[ll] = h__ * oldcs;
            e[ll]   = h__ * oldsn;
            if (*ncvt > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[nm13 + 1], &vt[ll + vt_dim1], ldvt);
            }
            if (*nru > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll * u_dim1 + 1], ldu);
            }
            if (*ncc > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[ll + c_dim1], ldc);
            }
            if ((r__1 = e[ll], std::abs(r__1)) <= thresh)
            {
                e[ll] = 0.f;
            }
        }
    }
    else
    {

        if (idir == 1)
        {
            f = ((r__1 = d__[ll], std::abs(r__1)) - shift)
                * (((d__[ll] > 0) ? c_b49 : -c_b49) + shift / d__[ll]);
            g    = e[ll];
            i__1 = m - 1;
            for (i__ = ll; i__ <= i__1; ++i__)
            {
                F77_FUNC(slartg, SLARTG)(&f, &g, &cosr, &sinr, &r__);
                if (i__ > ll)
                {
                    e[i__ - 1] = r__;
                }
                f            = cosr * d__[i__] + sinr * e[i__];
                e[i__]       = cosr * e[i__] - sinr * d__[i__];
                g            = sinr * d__[i__ + 1];
                d__[i__ + 1] = cosr * d__[i__ + 1];
                F77_FUNC(slartg, SLARTG)(&f, &g, &cosl, &sinl, &r__);
                d__[i__]     = r__;
                f            = cosl * e[i__] + sinl * d__[i__ + 1];
                d__[i__ + 1] = cosl * d__[i__ + 1] - sinl * e[i__];
                if (i__ < m - 1)
                {
                    g          = sinl * e[i__ + 1];
                    e[i__ + 1] = cosl * e[i__ + 1];
                }
                work[i__ - ll + 1]        = cosr;
                work[i__ - ll + 1 + nm1]  = sinr;
                work[i__ - ll + 1 + nm12] = cosl;
                work[i__ - ll + 1 + nm13] = sinl;
            }
            e[m - 1] = f;

            if (*ncvt > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "F", &i__1, ncvt, &work[1], &work[*n], &vt[ll + vt_dim1], ldvt);
            }
            if (*nru > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("R", "V", "F", nru, &i__1, &work[nm12 + 1], &work[nm13 + 1], &u[ll * u_dim1 + 1], ldu);
            }
            if (*ncc > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "F", &i__1, ncc, &work[nm12 + 1], &work[nm13 + 1], &c__[ll + c_dim1], ldc);
            }
            if ((r__1 = e[m - 1], std::abs(r__1)) <= thresh)
            {
                e[m - 1] = 0.f;
            }
        }
        else
        {

            f = ((r__1 = d__[m], std::abs(r__1)) - shift)
                * (((d__[m] > 0) ? c_b49 : -c_b49) + shift / d__[m]);
            g    = e[m - 1];
            i__1 = ll + 1;
            for (i__ = m; i__ >= i__1; --i__)
            {
                F77_FUNC(slartg, SLARTG)(&f, &g, &cosr, &sinr, &r__);
                if (i__ < m)
                {
                    e[i__] = r__;
                }
                f            = cosr * d__[i__] + sinr * e[i__ - 1];
                e[i__ - 1]   = cosr * e[i__ - 1] - sinr * d__[i__];
                g            = sinr * d__[i__ - 1];
                d__[i__ - 1] = cosr * d__[i__ - 1];
                F77_FUNC(slartg, SLARTG)(&f, &g, &cosl, &sinl, &r__);
                d__[i__]     = r__;
                f            = cosl * e[i__ - 1] + sinl * d__[i__ - 1];
                d__[i__ - 1] = cosl * d__[i__ - 1] - sinl * e[i__ - 1];
                if (i__ > ll + 1)
                {
                    g          = sinl * e[i__ - 2];
                    e[i__ - 2] = cosl * e[i__ - 2];
                }
                work[i__ - ll]        = cosr;
                work[i__ - ll + nm1]  = -sinr;
                work[i__ - ll + nm12] = cosl;
                work[i__ - ll + nm13] = -sinl;
            }
            e[ll] = f;

            if ((r__1 = e[ll], std::abs(r__1)) <= thresh)
            {
                e[ll] = 0.f;
            }
            if (*ncvt > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "B", &i__1, ncvt, &work[nm12 + 1], &work[nm13 + 1], &vt[ll + vt_dim1], ldvt);
            }
            if (*nru > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("R", "V", "B", nru, &i__1, &work[1], &work[*n], &u[ll * u_dim1 + 1], ldu);
            }
            if (*ncc > 0)
            {
                i__1 = m - ll + 1;
                F77_FUNC(slasr, SLASR)
                ("L", "V", "B", &i__1, ncc, &work[1], &work[*n], &c__[ll + c_dim1], ldc);
            }
        }
    }

    goto L60;

L160:
    i__1 = *n;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        if (d__[i__] < 0.f)
        {
            d__[i__] = -d__[i__];

            if (*ncvt > 0)
            {
                F77_FUNC(sscal, SSCAL)(ncvt, &c_b72, &vt[i__ + vt_dim1], ldvt);
            }
        }
    }

    i__1 = *n - 1;
    for (i__ = 1; i__ <= i__1; ++i__)
    {

        isub = 1;
        smin = d__[1];
        i__2 = *n + 1 - i__;
        for (j = 2; j <= i__2; ++j)
        {
            if (d__[j] <= smin)
            {
                isub = j;
                smin = d__[j];
            }
        }
        if (isub != *n + 1 - i__)
        {
            d__[isub]         = d__[*n + 1 - i__];
            d__[*n + 1 - i__] = smin;
            if (*ncvt > 0)
            {
                F77_FUNC(sswap, SSWAP)
                (ncvt, &vt[isub + vt_dim1], ldvt, &vt[*n + 1 - i__ + vt_dim1], ldvt);
            }
            if (*nru > 0)
            {
                F77_FUNC(sswap, SSWAP)
                (nru, &u[isub * u_dim1 + 1], &c__1, &u[(*n + 1 - i__) * u_dim1 + 1], &c__1);
            }
            if (*ncc > 0)
            {
                F77_FUNC(sswap, SSWAP)
                (ncc, &c__[isub + c_dim1], ldc, &c__[*n + 1 - i__ + c_dim1], ldc);
            }
        }
    }
    goto L220;

L200:
    *info = 0;
    i__1  = *n - 1;
    for (i__ = 1; i__ <= i__1; ++i__)
    {
        if (e[i__] != 0.f)
        {
            ++(*info);
        }
    }
L220:
    return;
}
